Issues and Challenges in Orbital-free Density Functional Calculations 
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Solving the Euler equation which corresponds to the energy minimum of a density functional 
expressed in orbital-free form involves related but distinct computational challenges. One is the 
choice between all-electron and pseudo-potential calculations and, if the latter, construction of the 
pseudo-potential. Another is the stability, speed, and accuracy of solution algorithms. Underlying 
both is the fundamental issue of satisfactory quality of the approximate functionals (kinetic energy 
and exchange-correlation). We address both computational issues and illustrate them by some 
comparative performance testing of our recently developed modified-conjoint generalized gradient 
approximation kinetic energy functionals. Comparisons are given for atoms, diatomic molecules, 
and some simple solids. 
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I. INTRODUCTION 

Investigation of orbital-free density functional theory 
(OF-DFT) 1-12 , including development of approximate 
orbital-free kinetic energy (OFKE) functionals, has at 
least two motivations. One is simply the beguiling no- 
tion of direct realization of the content of the Hohenberg- 
Kohn theoremi^r— . The other is practical, namely 
the possibility of eliminating the computational bottle- 
neck of solving the Kohn-Sham (KS) eigenvalue equa- 
tions, thereby dramatically broadening the applicability 
of Born-Oppcnheimer molecular dynamics run with DFT 
electronic energies. The practical aspect is the main fo- 
cus of the present work. 

In OF-DFT, the total electronic energy of an N e elec- 
tron system is a functional of the electron density n(r) 



£OF-DFTr , 



= T s [n] + #NeM + E H [n] 



E 



NN, 



(1) 



where T s [n] is the KS (non- interacting) kinetic energy 
functional given explicitly as a density functional, E^ c [n] 
is the nuclear-electron interaction energy, .EhM is the 
Hartree energy (classical electron-electron repulsion), 
E xc [n] is the exchange-correlation (XC) energy func- 
tional, and -Enn is the inter-nuclear repulsion energy. 
Minimization of the functional Eq. ([1]) gives a single Eu- 
ler equation to be solved, 



ST s [n] 
5n(r) 



+ «Ks(M; r ) = m- 



(2) 



Here vks is the Kohn-Sham potential, S(E-^ C + Ejj + 
E xc )/Sn and /j, is the chemical potential. The ordinary 
KS equation has the same potential but requires solu- 
tion for N e or N e /2 orbitals (in the all-electron, spin- 
polarized and non-spin-polarized cases respectively) . So- 
lution of the ordinary KS problem scales computationally 
as w iVg in general, whereas solution of Eq. ([5} should 
scale approximately linearly. 



Practical implementation of OF-DFT requires approx- 
imation of both T s [n] and E xc [n] . Simply because of their 
relative magnitudes, the quality of an OF-DFT calcula- 
tion is dominated by the quality of the approximate T s . 
There are two distinct classes of approximation in the 
literature, one-point functionals, 



T s [n] = / t s ([n];v)d 3 



(3) 



and two-point functionals 

T s [n] = // 1)S ([n];r)x(r ) r')/ 2 , s ([n];r')d 3 rd 3 r'. (4) 

Here /i jS and /2, s are weighting functionals and x( r , r ') 
is a type of response function. For reasons of compu- 
tational efficiency as well as conceptual simplicity (two- 
point functionals take the development out of the frame- 
work of an effective Kohn-Sham equation (see Eq. ([8|) 
below) unless an optimized effective potential^ is used, 
itself an extra complication), we (and our collaborators) 
have focused exclusively on one-point functionals and do 
so here as well. 

An interesting feature of the literature on developing 
approximate OFKE functionals, including our contribu- 
tions with collaborators, is that there are more tests of 
approximate functionals using inputs from other sources 
(e.g. conventional KS calculations, Hartree-Fock calcu- 
lations, etc.) than tests by solving the Euler equation, 
Eq. ©. A side effect is that comparatively little is known 
about the difficulty of solving that equation with approxi- 
mations other than of the Thomas- Fermi kind (see below) 
and about the relative effectiveness of various solution 
techniques. 

To frame that issue and the calculations reported here, 
it is useful to decompose the non-interacting KE func- 
tional into the von Weizsacker contribution^ plus a non- 
negative remainder, the Pauli termor— , 



T B [n]=T w [n]+T e [n], T e [n] > 



(5) 



2 



The von Weizsacker functional (in Hartree atomic units) 



Twin] = - 



1 f \Vn(r)\< 



8 J n(r) 



d 3 r = / t w {[n};r)d 3 r . (6) 



It is exact for one electron and for a two-electron singlet. 
From 



ST w [n] 1 . 1 



Sn(r) 2 



(--V 2 )^W, 



(7) 



the Euler equation Eq. ([TJ) takes a Schrodinger-like 
formS 3 .^ 

1 



V 2 +v e ({n];r) +v KS ([n};r) \ y/n(r) = MV^W 



(8) 

Observe that, unlike familiar quantum mechanical eigen- 
value problems, the "orbital" in Eq. (j8|) is normalized to 
N e , not unity. Here vg is the Pauli potential, 



v e ([;;]; r) 



ST s [n] 



Sn(r) 



vg([n};r) >0 . 



(9) 



Non-negativity of Tg and vg has proved to be an 
important pair of constraints for OFKE functional 
development^— . 

Eq. ([8]) resembles the ordinary KS equation, a fact 
that has led to contradictory statements about solution 
techniques. On the one hand, Ref. declares that Eq. 
(|5]) "... can be solved iteratively to self-consistency by 
any Kohn-Sham computer program: just select the low- 
est eigenvalue. The solution is very simple and quick, 
for there is only one 'orbital' . . . " . Ref. [26| makes 
precisely the contrary claim, at least in the context 
of the widely used Gaussian-type orbital (GTO) basis 
sets. Those authors expanded ^/n in a GTO basis with 
coefficients Cj, with respect to which they minimized 
C := E OF - DFT [n] - nN e . They state that "Due to the 
highly nonquadratic nature of the kinetic energy, the op- 
timization of £ is a nontrivial problem. The iterative 
self-consistent procedure used in Kohn-Sham calculations 
does not work, and we require more robust minimization 
techniques. Moreover, . . . first derivative methods such 
as conjugate gradient minimization and quasi-Newton 
search perform poorly, requiring many hundreds of it- 
erations to achieve convergence." A related discussion 
and references to the few earlier papers on the issue is at 
p. 135 of Ref. [TJ. This is one of the issues addressed in 
the present study. 



II. APPROXIMATE KINETIC ENERGY 
FUNCTION ALS 

To set the stage for another technical issue, we consider 
types of approximate one-point OFKE functionals next. 



For work on minimization involving two-point function- 
als, see Refs. 31, 32 and references therein. 



A. Thomas-Fermi Type 

Diverse approximate OFKE functionals can be written 
in the generic form 

T s [n] =Tw[n] + ATtfM + T A [n] 

< A < 1 . (10) 

The simplest local approximation for the KE is the 
Thomas-Fermi (TF) 33 ' 34 functional 

T TF [n}= J t TF ([n};r)d 3 r = c J n 5 / 3 (r) d 3 



Co 



(ii) 

alone. The approximation Ta = and A = 1 is widely 
used in many OF-DFT applications (see Ref. for discus- 
sion and references) despite its known deficiencies^. A 
related form, commonly called Thomas-Fermi-Dirac-von 
Weizsacker theory, is a linear combination of Ttf with 
some fraction of Tw, 



TFvW.q 



w 



■T- 



TF 



< a < 1 



(12) 



along with the local Dirac exchange functional. Early re- 
ports of special self-consistent OF-DFT calculations men- 
tioned above were for this mode l 26 ' 36 ' 37 . 

As an aside, there is an extensive literature of ef- 
forts to determine an optimal value of a in Eq. (|12jl . 
Since the Pauli term decomposition, Eq. ([5]), provides 
both an exact lower bound on T s and leads to the den- 
sity Schrodinger equation, Eq. ([8]), that decomposition, 
and its elaboration Eq. (fTt)|) . seems preferable to using 
TtfvW,q; and attempting to optimize a. But because 
TtfvW,o: is prevalent in the literature, we consider the 
numerical issues associated with it as well. 



B. Generalized Gradient Approximation KE 
Functionals 

Generalized gradient approximations (GGA) are best 
known in DFT as improvements on the local approxima- 
tion for E xc . For either E xc or T s , a GGA is a truncation 
of the corresponding gradient expansion which is altered 
to meet relevant constraints and suppress unphysical be- 
haviors. For the KE functional, a GGA can be written 
as 



?f GA M= / t TF ([n];r)F t ( S (r))d 3 r , 



(13) 



where Ft is the kinetic energy enhancement factor. It is 
a function of the dimensionless reduced density gradient, 

\Vn\ I |Vn| (14) 



(2k F )n 2(3tt 2 ) 1 /3 n 4 / 3 
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Because f\y 



|s 2 t T F, the GGA Pauli term in Eq. (0 is 



T, 



GGA r 



*] = / tT F ([n];r)F e (s(r))d 3 



(15) 



Ref. llCj showed that the KS KE of a molecular sys- 
tem is dominated by the behavior of Fg over a relatively 
small range of s. For much of that range, Figure [T] dis- 
plays the Pauli enhancement factors for the functionals 
TtfvW,q, with a = 1, 1/9 Eq. ([12]). the Tran-Wesolowski 
GGA22 .(PBE-TW), and the mcGGA functional (PBE2) 
of Ref. S The latter two use the same enhancement fac- 
tor form as the Perdew, Burke, and Ernzerhof (PBE) 39 
GGA X functional, F x (s) = 1 + cs 2 /(l + as 2 ). In PBE- 
TW F t oc F Xj pbe with parameters fitted to reproduce 
the kinetic energy of a small training set, an assumption 
called conjointness. PBE2 is a "modified conjoint" GGA 
(mcGGA) functional because the parameters in it were 
constrained to satisfy Pauli-term non-negativity; see Ref. 
M for details. 

Observe in Fig. [TJ that the PBE-TW and T TFv w !Q=1/9 
Pauli enhancement factors are almost identical, espe- 
cially for small s. There, both have negative slope 
(with respect to s 2 ) which causes violation of vg non- 
negativity^, recall Eq. ([§]). The common property of 
the TrFvW,a=i an d PBE2 approximations is satisfac- 
tion of that non-negativity constraint. The low slope 
of the PBE2 enhancement factor at small values of s 2 , 
F e FBE2 (s) w 1 + 0. 3642s 2 , makes the enhancement factors 
for Tt F vw,q=i and Tpbe2 close for s < 1. This compari- 
son suggests that the results obtained with the PBE-TW 
KE functional should be close to those from 7T FvWQ , =1 /g 
and, similarly, the results from PBE2 should be close to 
those from TrFvW,a=i- 

A technical problem common to these GGAs is that 
both fe, pbe-tw and vg,pBE2 are singular at nuclear sites, 
the former negative, the latter positive. Numerical so- 
lution of the Euler equation Eq. |2]) must address this 
problem, an issue to which we return below. 



III. ALL-ELECTRON SOLUTIONS OF THE 
OF-DFT EULER EQUATION 

As in ordinary KS calculations, solution of Eq. ([5]) can 
be either all-electron or via pseudo-potentials. In this 
section, we consider all-electron solutions, by both GTO- 
basis and numerical grid techniques and address pseudo- 
potentials in the subsequent section. 



A. Atoms 

To test the notion that any standard KS code can 
be used straightforwardly 24 , we modified the GTO-basis 
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FIG. 1: Pauli term enhancement factors Fg of OFKE func- 
tionals as a function of s 2 . GGA denotes the Tran-Wesolowski 
functional, mcGGA denotes the PBE2 functional. See text for 
details. 



code SOAtom to handle Tg and vg as in Eqs. (0, ©, 
and ©. SOAtom, a part of the GTOFF suite^MI, solves 
the KS equation in a Hermite Gaussian basis with ana- 
lytical evaluation of all the matrix elements except for 
those involving XC. Those are done on a radial grid. 
We also modified the Laaksonen all-numerical diatomic 
molecular code^ 2 - correspondingly. It is based on a pro- 
late spheroidal grid. 

Insofar as numerical stability is concerned, the results 
are quite clear. Even for TxFv\v,a=i with simple Slater 
exchange (i.e. TFvWD), the typical iterative SCF pro- 
cedure is only marginally stable. The problem is worse 
in the GTO basis than in the grid-based calculation, at 
least in the specific sense that a simple SCF stabiliza- 
tion procedure (Pratt, i.e. linear mixing of a fraction 
of current iteration density and the rest from the preced- 
ing iteration) fails completely for many OF-DFT calcula- 
tions. Ordinary KS calculations on the same atoms with 
the same simple stabilization scheme converge in a few 
iterations. 

The lithium and carbon atoms are examples. For Li in 
a 9s GTO basis in the SOAtom code, the pure TF form 
(T B = T TF , i.e., Eq. (TJJ) with a = 0), the T e = T TF 
form (Eq. (TO} with A = 1, T A = 0), and the Tg = AT TF 
form (Eq. (flTJ|) with Ta = 0) can be brought to numerical 
convergence but the mcGGA form Tg = T mc QGA — T v w 
cannot. For the successes, more iterations by one to 
two orders of magnitude are required than for conven- 
tional KS and the numerical convergence is poor. One 
can get to fractional total energy errors of 10~ 4 — > 10~ 2 
for lighter to heavier atoms respectively, but not much 
better. The contrast with conventional KS atomic calcu- 
lations is stark: in them convergence to 10 -6 is trivial to 
achieve. 

Table U illustrates this point with comparison of nu- 
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merical grid and GTO-basis results for T s = T v \v and 
T v w + Ttf with simple Slater exchange (Dirac exchange) 
on H, Li, and Ne. (The GTO calculations are with 9s 
basis sets for H and Li, 13s for Ne.) The two total ener- 
gies for Li differ at the 1 mHartree scale. Notice that the 
numerical-grid results match rather well with the values 
from Ref. l26l which were calculated with a direct min- 
imization scheme, not a modified KS code. Ironically, 
a misbehavior of simple Slater exchange, namely that it 
satisfies the virial theorem in the form E to t = —T s (which 
the exact E xc does not), in this case highlights the con- 
vergence problem, especially in the GTO calculation. 

Results for the carbon atom in the GTO basis, not 
shown in the Table, are worse. The TTFvW,a=i/5 calcu- 
lation with a 13s basis cannot be brought to SCF conver- 
gence, even with tricks such as starting with full T\y and 
no Ttf contribution, then slowly scaling down the former 
while scaling up the latter. The corresponding standard 
KS calculation converges trivially. 

For the numerical-grid calculations, SCF convergence 
is very slow compared to standard KS calculations, but 
reasonable results can be obtained. Table |TT] shows total 
energies for the first row atoms obtained from numerical- 
grid self-consistent OF-DFT calculations with various 
OFKE approximations, again with Slater exchange. For 
TFvWD and TrFvW.r>= i /b l /g, comparison with the di- 
rect minimization of Ref . |2£j| (the first six columns of data) 
confirms that our calculations succeeded. 

Note that the total energies from the TF+vW and 
mcGGA(PBE2) kinetic energy functionals are overesti- 
mated (as a consequence of overestimation of the KS 
KE). In contrast, all of the functionals with scaled von 
Weizsacker contributions underestimate the KS KE, so 
that the resulting total energies are below the reference 
KS values. Such behavior is characteristic of a failure of 
iV-representability in the KE functional^ 3 -. Observe also 
that TF+l/9vW and GGA(PBE-TW) total energies are 
close to each other, though the functional forms differ. 

Table HTT1 shows the effects of using the full LDA E xc , in 
this case the VWN parameterization^. Unsurprisingly 
but reassuringly, inclusion of the C functional shifts the 
total energies downward without altering the trends. 

These atomic results lead us to nuanced agreement 
with the claim of Ref. [26| and disagreement with the claim 
of Ref. 24. The OF-DFT Euler equation is not, in gen- 
eral, solvable by simple modification of a standard GTO 
KS code (the norm for molecular calculations). Even a 
good all-numerical KS code is challenged to achieve so- 
lutions but can be made to succeed for isolated atoms. 
Realizing the computational speed-up potential of OF- 
DFT clearly depends on algorithms and implementations 
well-suited for OF-DFT, even for one-point functionals. 



B. Diatomic Molecules 

Numerical-grid solution of Eq. © for diatomic 
molecules, if possible, would yield two kinds of insight: 
numerical method behavior and the comparative behav- 
ior of n(r) and vg(r) for different OFKE approximations. 
Though the difficulties of using a modified KS code are 
just as evident in this case, we have been able to achieve 
solutions for several light molecules. 

Numerical requirements include extremely tight con- 
vergence tolerances on the eigenvalue fi and normaliza- 
tion (10~ 3 more stringent than normal KS calculations), 
much larger maximum distance cutoff (80 to 100 a.u. 
vs. normal KS 30 to 40 a.u.), and about a factor of five 
more points in both of the prolate spheroidal coordinates 
(roughly 1100 x 1300 points vs. the typical 200 x 300). 
Even so, the total energy convergence is mediocre for 
GGA and mcGGA functionals, about 0.01 Hartree at 
best. Convergence is better for the TtfvW.q function- 
als, between 0.1 and 1 mHartree. This need for extreme 
measures to achieve limited-quality outcomes is an ad- 
ditional confirmation of the unsuitability of unmodified 
conventional KS schemes for solutions of the OF-DFT 
Euler equation. 

The solutions nevertheless provide real compar- 
ative insight regarding different OFKE approxima- 
tions. Figure [5] compares the all-electron KS density 
(£"xc,lda, VWN) around the Si site in SiO with the 
densities from Tgga,pbe-tw (the Tran-Wesolowski 38 
GGA), T mcGGA ,PBE2 (the PBE2 mcGGA 9 ), and 
TtfvW,q=i/9,i- These approximate functionals form 
pairs. TggAjPBE— tw pairs with r T Pvw,a=i/9) while 
TtFvW,c*=i pairs with T mc GGA,PBE2- This pairing con- 
forms to the expectations formed in considering the 
small-s behavior of the respective enhancement factors. 
The pairing also is interpretable directly from the near- 
nucleus repulsion or attraction behavior of the vari- 
ous approximations. Tgga,pbe-tw generates a vg with 
a spurious negative singularity near the nuclei, while 
^TFvW,o!=i/9 drastically lowers the von Weizsacker lower 
bound to T s . Both lead to excess near- nucleus density. In 
contrast, we, m cGGA,PBE2 has spurious positive nuclear site 
singularities^. Near the nuclei, however, f e,mcGGA,PBE2 
and i>0,tfvW,q=i match quite well, as shown in Fig. [3] 
The result, shown in Fig. [51 is that these two functionals 
give rather close to the same density. Observe that the 
behavior of vg in the vicinity of the nucleus for each of 
these approximate functionals differs dramatically from 
that of vg obtained by inversion of the standard KS 
scheme. It is that improper behavior which we believe 
causes the problems with convergence of standard KS 
codes used with approximate OFKE functionals. The 
positive near-nuclei singularities appear, in particular, to 
pose numerical problems. 

Details of the density near the Si site are provided 
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TABLE I: All-numerical and GTO results for the atoms H, Li, and Ne for the TrFvW,a=o,i models with simple Slater exchange. 
Energies in Hartree a.u. 

T s = T w Numer. T s = T w GTO T s = T w + T TF Numer. T s = T w + T TF GTO T s = T w + T TF a 



H Atom 

E ta t -0.406534 -0.400737 -0.261827 -0.259969 -0.2618 

T s 0.406534 0.859699 0.261827 0.262042 

T g 0.000 0.000 0.091034 0.090221 

ji -0.1943 -0.1764 -0.0715 -0.0696 -0.071 

Li Atom 

Etot -8.525825 -8.523413 -4.105425 -4.096347 -4.1054 

T s 8.525825 8.523126 4.105425 4.103660 

T g 0.000 0.000 2.019249 2.009622 

£ -0.9575 -0.9526 -0.1306 -0.0.1365 -0.131 

Nc Atom 

Etot -274.68080 -274.652253 -85.734451 -85.730041 -85.7343 

T s 274.68080 274.664688 85.734438 85.728273 

Tg 0.000 0.000 54.352106 54.347495 

u -7.0607 -7.0594 -0.1807 -0.1806 -0.181 



"From Ref. [26| 



TABLE II: Self-consistent atomic total energies obtained from various OFKE functionals (Hartree a.u.) and simple Slater 
exchange. 



1/9 vW+TF" 1/9 vW+TF 1/5 vW+TF" 1/5 vW+TF vW+TF a vW+TF GGA (PBE-TW) mcGGA (PBE2) KS b 



H 


-0.6664 


-0.6664 


-0.5666 


-0.5666 


-0.2618 


-0.2618 


-0.71 


-0.32 


-0.4065 


He 


-3.2228 


-3.2228 


-2.8184 


-2.8184 


-1.4775 


-1.4775 


-3.4 


-1.5 


-2.7236 


Li 


-8.2515 


-8.2515 


-7.3227 


-7.3227 


-4.1054 


-4.1054 


-8.6 


-4.1 


-7.1749 


Be 


-16.1631 


-16.1631 


-14.4841 


-14.4841 


-8.4922 


-8.4922 


-16.7 


-8.4 


-14.2233 


B 


-27.2876 


-27.2876 


-24.6284 


-24.6284 


-14.9258 


-14.9259 


-28.0 


-14.6 


-24.5275 


C 


-41.9052 


-41.9053 


-38.0332 


-38.0332 


-23.6568 


-23.6569 


-42.9 


-23.0 


-37.6863 


N 


-60.2622 


-60.2623 


-54.9428 


-54.9429 


-34.9084 


-34.9084 


-61.6 


-33.9 


-54.3977 


O 


-82.5798 


-82.5799 


-75.5765 


-75.5765 


-48.8831 


-48.8832 


-84.3 


-47.3 


-74.8076 


F 


-109.0592 


-109.0594 


-100.1345 


100.1346 


-65.7674 


-65.7676 


-111.1 


-63.5 


-99.4072 


Ne 


-139.8865 


-139.8867 


-128.8014 


-128.8016 


-85.7343 


-85.7344 


-142.3 


-82.7 


-127.4907 



"From Ref. |26|. 

'Spin-restricted LDA (Slater exchange) calculation. 



in Fig. |U For purposes of display, the densities are 
weighted by a quasi-radial factor with origin at the Si 
site, 47r(|z| — R/2) 2 . The proper KS shell structure is 
missing, as is usual with single-point OFKE functionals. 
The more repulsive nature of the pair T' mc GGA,PBE2 and 
TtfvW.q^i compared to T G ga,pbe-tw and T TFvW ,a=i/9 
also is evident. It is interesting that in the region 
—2.5 < z < 2.0 T mc GGA,PBE2 does give a weak mimicry 
of the outermost shell structure in T s , unlike the other 
models. We are uncertain as to how reliable or useful 
this feature is. 

Fig. [5] compares the behavior of E OF ~ BFT [n] as a func- 
tion of SiO bond length with standard KS results. One 
sees immediately that the GGA and mcGGA forms intro- 
duce numerical difficulties because of their dependence 
on the reduced density gradient s, Eq. (p~4|) . Clearly 
there is a grid interval-size problem which could be obvi- 



ated by going to even denser grids but at obvious com- 
putational cost. The known failure of Tgga.pbe-tw to 
give binding^ is evident. i?TFvW,a=i/9 apparently does 
not bind either, in keeping with the too- weakened lower 
bound just discussed. Full TFvWD and £"mcGGA,PBE2 
are fairly close, with the mcGGA being the best of the 
lot with respect to equilibrium bond length. 



C. Simple Analysis of the Difficulty 

The barrier to use of a standard KS code to solve Eq. 
([5]) can be traced to the near- nucleus repulsion of vg . As 
displayed in Fig. [31 the exact vg is strongly repulsive in a 
fairly small region around the nuclear site. Ref. El Fig. 
2 shows that the exact vg can have rather sharp struc- 
ture within a radius of about 1 bohr of a nuclear site. 
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TABLE III: OF-DFT self-consistent atomic total energies (Hartree a.u.) obtained from various kinetic energy functionals and 
VWN -E/xc.lda with the numerical grid KS code. 





1/9 vW+TF 


1/5 vW+TF 


vW+TF 


GGA (PBE-TW) 


mcGGA (PBE2) 


KS a 


H 


-0.7101 


-0.6084 


-0.2924 


-0.76 


-0.36 


-0.4457 


He 


-3.3244 


-2.9175 


-1.5590 


-3.5 


-1.6 


-2.8348 


Li 


-8.4175 


-7.4860 


-4.2469 


-8.7 


-4.2 


-7.3352 


Be 


-16.3982 


-14.7162 


-8.6995 


-16.9 


-8.5 


-14.4472 


B 


-27.5953 


-24.9329 


-15.2033 


-28.4 


-14.8 


-24.3436 


C 


-42.2886 


-38.4132 


-24.0078 


-43.3 


-23.4 


-37.4202 


N 


-60.7237 


-55.4007 


-35.3357 


-62.1 


-34.3 


-54.0250 





-83.1215 


-76.1146 


-49.3893 


-84.8 


-47.8 


-74.4613 


F 


-109.6832 


-100.7547 


-66.3545 


-111.7 


-64.1 


-99.0960 


Ne 


-140.5945 


-129.5054 


-86.4042 


-143.1 


-83.4 


-128.2335 



"Spin-restricted LDA calculation. 
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FIG. 2: All-electron self-consistent Kohn-Sham and OF-DFT 
electron densities plotted along the SiO molecule internuclear 
axis in the vicinity of the Si site. Si at (0,0,-1.05) A, O out of 
the picture at (0,0, +1.05) A. See text. 
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FIG. 3: Pauli potentials ve around the Si site in SiO from self- 
consistent all-electron Kohn-Sham and OF-DFT calculations 
with the PBE2 mcGGA, TW GGA, and T TFvWia=1/9il OFKE 
functionals. 




-3 -2.5 -2 -1.5 



z (angstrom) 

FIG. 4: All-electron self-consistent Kohn-Sham and OF-DFT 
electron densities near the Si site along the SiO molecular 
axis. These are scaled scaled by the factor 47r(|z| — R/2) , 
with R = 2.10 A, the internuclear distance. This puts the 
origin of the scaling at the Si site (0,0,-1.05) A. The O is out 
of the picture at (0,0, +1.05) A. 



In contrast, some simple approximations which are prop- 
erly positive definite, including our mcGGA, actually are 
singular at the nuclei; again see Fig. [3l Such strong re- 
pulsion overwhelms the attractive v xc . That figure also 
shows that some approximations deliver Pauli potentials 
with negative nuclear-site singularities. We consider that 
case below. First, however, the simplest example will suf- 
fice to illustrate the problem with properly positive vg. 
Pick Tg = Ttf, Eq. (fTTj). and E xc to be simplest Slater 
exchange: 

«=-!(}) " 3 (16) 
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1 1.5 2 2.5 

R (angstrom) 

FIG. 5: Total energy of the SiO molecule as a function of 
bond length obtained from self-consistent all-electron Kohn- 
Sham and OF-DFT calculations with Thomas-Fermi, Tran- 
Wesolowski (GGA) and PBE2 (mcGGA) kinetic energy func- 
tionals. Kohn-Sham values are shown for comparison. Val- 
ues are shifted to a common zero by 363.076 (KS), 386.339 
(TF+vW/9), 250.529 (TF+vW), 83.902 (GGA) and 241.564 
(mcGGA) Hartree a.u. 



Then in Eq. ((5]), the potentials become 
v e =i(3^ 2 ) 2 /V/ 3 (r) 



1/3 



A hydrogen-like density, 



n H (r) 



exp(-2A^ e r) 



(17) 



(18) 



obeys the Kato cusp condition near the nucleu o 16 ' 27 " — , 
hence is useful for testing. At the nucleus, this density 
yields the ratio of potentials 



Pxc(0)+P 9 (0) 

K(o)l 



= -l + 3.318004iV p 4/3 



(19) 



For N e — 6, this ratio is already 35.2. By N e — 10 it is 
70.5. Additional simple calculations with the potential 
which appears in Eq. ((5J without the positive Hartree 
contribution, that is VN e (r) +v$(r) + v KC (r), illustrate the 
point. At small r with N e = 6, that potential becomes 
positive for r > 0.028 Bohr. For N e = 10 it is positive 
for r > 0.011 Bohr. These little exercises illustrate why 
the use of an ordinary K-S code becomes so difficult. 
Such peculiar behavior is quite different from what is 
encountered in the ordinary KS problem. 

From the perspective of numerical stability, the case 
of a vg which has a spurious negative singularity near 



each nuclear site, e.g. PBE-TW, is at least as bad if 
not worse. As a site is approached, such potentials first 
are increasingly repulsive, then plunge abruptly into the 
negative singularity; see Fig. [3] 



IV. PSEUDO-POTENTIAL SOLUTIONS OF 
THE OF-DFT EULER EQUATION 

Having demonstrated the difficulties with solving the 
OF-DFT Euler-Lagrange problem with a modified KS 
eigenvalue code, we turn to the use of direct Euler- 
Lagrange minimization. The specific objective is to 
exploit the numerical methodology in the Profess 
cod o 31 ' 32 . Written originally for use with two-point func- 
tionals, Profess performs OF-DFT calculations by min- 
imization of the Euler-Lagrange equation as a functional 
of n(r) under periodic boundary conditions. It uses a 
numerical 3D mesh and FFTs. As published, the code 
includes the TF, vW, and TFvW, a functionals as well 
as the Wang-Teter (WT)^, and Wang-Govind-Carter 
(WGC)^ OFKE functionals. The PZ and PBE E xc func- 
tionals are implemented in Profess. For this study, 
we added the Tran-Wesolowski GGA 3 ^ and our PBE2 
mcGGA 9 OFKE functionals. 

As is the case with standard KS calculations done in a 
plane-wave basis, Profess relies upon pseudo-potential 
(PP) techniques to screen the nuclear-electron potential 
cusp and exclude chemically inactive core states. Though 
OF-DFT has no problem with core states and the density 
(and its square root) is a comparatively unstructured, 
smooth function, regularization of the nuclear-electron 
interaction singularity still is a requisite for an efficient 
implementation. 

High-quality pseudo-potentials developed for conven- 
tional KS calculations generally are non-local, in the spe- 
cific sense that they contain projection operators which 
provide different potentials for different orbital angular 
momenta. That explicit orbital dependence makes non- 
local pseudo-potentials (NLPP) inapplicable in OF-DFT 
calculations. Instead, local pseudo-potentials (LPP), i.e., 
of the form of a simple multiplicative operator which is 
the same for all orbitals, must be developed. Profess 
requires a LPP in real or reciprocal space as input. Ob- 
serve that this limitation to local form is an additional 
approximation, over and beyond the PP itself, which has 
accuracy limitation implications for both conventional, 
orbital-based KS or OF-DFT implementations. 

In addition to their simplicity, there is a formal ad- 
vantage of LPPs which is at least of peripheral interest 
here. Calculations with local PPs are within the frame- 
work of the standard KS scheme, which assumes a local 
effective potential. The NLPP case obviously does not 
meet that assumption. Although the Hohenberg-Kohn 
theorem has been extended to the case of a non-local 
external potential 4 ^, the exchange-correlation energy in 
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that case becomes be a functional of the one-particle re- 
duced density matrix instead of a functional n(r) alone. 

Many methods have been proposed to develop LPPs. 
Among them we mention (i) empirical (or model) LPPs 
as, for example in Refs. I48l452t (ii) local potentials ob- 
tained from non-local ones, for example, by use of just 
one Z-channel from an NLPP as a LPP (for example, 
Ref. [53|; (iii) LPPs constrained to reproduce atomic prop- 
erties, eigenvalues, or pseudo-density, etc., which follow 
from a (presumably superior) NLPP (for example, Ref. 
54[ ) , and finally (iv) local PPs derived to reproduce some 
bulk property values, either experimental or those pre- 
dicted by NLPP calculations^^ 



A. Local Pseudo-potentials for OF-DFT 
Calculations 

1. Development 

Some time ago, an iterative procedure was 
develope d 57 ' 58 to solve the inverse problem of de- 
termining the KS effective potential t>Ks( r ) from a given 
density n(r). Subsequently, we^ introduced and tested 
an improvement. In the case of Li, however, both ver- 
sions share a problem. For a single valence orbital (singly 
or doubly occupied) the solution of the inverse problem 
is trivial and known. The local pseudo-potential is 
equal to the s-channel of the NLPP, vi oca i(r) — u; = o(r). 
Hence the LPP contains no information about the I > 
channels of the NLPP. Those channels are critical in 
crystalline binding. 

Therefore, to include information about all I channels 
of the reference NLPP, we consider a sort of normalized 
linear combination of I components of that NLPP, 



,(r) = civi(r) Y 



ci 



(20) 



1=0 



1=0 



where the parameters {c;} are to be adjusted to fit se- 
lected equilibrium bulk material properties calculated 
with the reference KS method. This particular method 
of LPP generation amounts to a mixture of methods (ii) 
and (iv) described at the outset of this Section. 

In the present case, we simply took the bec Li lattice 
constant as predicted by a standard KS calculation with 
PBE 39 E*c, and the plane wave (PW) basis set (see Table 
TV|) . namely a — 3.44 A. Components of the Troullier- 
Martins norm-conserving NLPP were used in Eq. ([2"Tj)) . 
For generation of the NLPP with PBE XC, we took the 
core radius to be 2.45 a.u. The parameters {c/} in Eq. 
(|20p . for the s, p, and d channels respectively, were deter- 
mined by constraining a KS calculation with the LPP Eq. 
(|20| to reproduce the reference optimized bec Li lattice 
constant value. Those KS calculations done with PBE 
XC in the Siesta code 62 and a DZP numerical atomic 



orbital (NAO) basis set. The optimized parameter val- 
ues are cq = 0.69, ci = 0.34, C2 = 0.10. We designate 
this LPP as «GGA,spdi- To generate the LDA local LPP, 
WLDA, sp di for the Perdew-Zunger— LDA XC functional, 
components of the LDA NLPP and the same set of the 
channel-mixing parameters were used in Eq. (|20j) . 

An alternative LPP form which we also studied is 
a modification of the potential proposed by Heine and 
Abarenko v 49 ' 51 . In real space, the Heine- Abarenkov 
model potential is 



Vmod(r) 



-A, r < r c 
-Z/r, r >r c 



(21) 



where A is a constant, r c is the core radius, and Z is the 
core charge. The model potential in reciprocal space is 
given by 



^mod 



-4-7T 

TV 



[(Z-Ar c )cos(qr c ) + (A/q)sm(qr c )}, (22) 



where O is the unit cell volume. In Ref. |5lJ, this 
potential was multiplied by a smoothed step function 
/(<?) = ex P[ — q/Qc) 6 } to reduce spurious oscillations in 
v mo d{q) and to ensure rapid decay of v mo <i(q) at large 
wave- vectors. Those oscillations are caused by the dis- 
continuity of the real-space potential at the core radius. 
Here, the parameter q c was chosen as suggested in Ref. 
5ll . namely, to equal the second zero position of i> mo d(<z)- 

To obtain counterparts of the local potentials 
«GGA,s P di and «LDA,spdi, Eq. in the simple mod- 
ified Heine-Abarenkov model form, the two parame- 
ters, A and r c , were determined by minimization of 
/ dr|v S pdi(r) - v mo di(r)| 2 - This yields A = 0.45499 
Hartrees, r c — 2.2261 Bohr, g c = 2.86 Bohr -1 for 
fGGA.modi and A = 0.45376 Hartrees, r c = 1.8818 Bohr, 
g c = 2.94 Bohr -1 for f LDA,modi ■ The local potentials 
WGGA.spdi and i>LDA,spdi in reciprocal space are multi- 
plied by the same smoothed step function f(q) with q c 
values equal to 2.95 and 3.25 Bohr -1 respectively. 

Figure [5] shows the WGGA.spdi LPP in real space in 
comparison with the NLPP I channels, along with the 
two pseudo-densities which result. Figure [7] shows the 
«GGA, sp di and i>GGA,modi LPPs in reciprocal space. 



2. KS Tests of Local Pseudo-potentials 

Kohn-Sham calculations with the LPPs were per- 
formed using the Abinit PW code 64 with PZ and PBE 
exchange-correlation functionals. We also used Siesta 
with the same exchange-correlation functionals and a 
2s 2 2p 2 numerical atomic orbital basis set (8 NAO per 
atom). Table HVl shows the equilibrium lattice constants 
and bulk moduli for the various LPPs. Those results are 
compared to the Kohn-Sham calculations performed with 
the non-local projector augmented wave (PAW) scheme 
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FIG. 6: Real space pseudo-potentials for Li: local UGGA, sp di, 
and different /-components of the non-local Troullier-Martin 
(TM) pseudo-potential. Pseudo-densities generated with lo- 
cal and non-local PPs are shown for comparison. 
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FIG. 7: Reciprocal space local pseudo-potentials for Li: 

«GGA,spdl and «GGA,modl- 



(as implemented in Vasp and Abinit) and TM norm- 
conserving pseudo-potentials with core correction^. The 
lattice constant and bulk modulus reported in Table HVl 
were obtained by fitting the calculated total energies 
per cell to the stabilized jellium model equation of state 
(SJEOS^). All the local PPs reproduce the PAW results 
rather closely for both lattice constant and bulk modu- 
lus. The bulk moduli calculated using NAO orbitals and 
norm conserving TM pseudo-potentials are slightly larger 
than the PAW plane wave results. 

As a check against an all-electron localized-orbital cal- 
culation, we did high-quality GTO-basis KS calculations 
(10s6p3d basis) with the GTOFF code 41 . For E XC , PZ and 
-Exc,pbe, we obtained optimized bcc Li lattice parameters 
of 3.360 and 3.435 A, respectively, essentially the same as 
from the Siesta NAO and plane wave PAW calculations. 



TABLE IV: Kohn-Sham lattice constant (A) and bulk mod- 
ulus (GPa) for bcc Li calculated using Vasp PW PAW 
schemes, Abinit PW PAW and local pseudo-potentials, 
Siesta non-local Troullier-Martins^ and local pseudo- 
potentials. Orbital-free calculations used T m cGGA,PBE2, 
TtfvW,<*=i, Tgga,pbe-tw, and T TPvW>Q =i /9 kinetic energy 
functionals in combination with -E xc ,lda,pz and £ X c,gga,pbe 
with local pseudo-potentials «LDA,s P dl, «GGA,spdi, fLDA,modi 
and WGGA.modi- Conventional KS calculations were done with 
a 2-atom unit cell and 7x7x7 (Vasp and Siesta) or 9 x 
9x9 (Abinit) k-mesh. The Siesta basis set was 2s 2 2p 2 (8 
NAO per atom). Orbital- free calculations used a 128-atom 
supercell. 



LDA GGA 



Method 


PP 


a 


B 


a 


B 


Kohn-Sham 












PW (Vasp) 


PAW 


3.37 


15.0 


3.45 


13.7 


PW (Abinit) 


PAW 


3.37 


15.1 


3.44 


13.9 


NAO (Siesta) 


TM 


3.37 


15.6 


3.44 


14.3 


Kohn-Sham 












PW (Abinit) 


spdl" 


3.37 


14.8 


3.44 


13.8 


NAO (Siesta) 


spdl" 


3.38 


14.9 


3.45 


13.9 


Kohn-Sham 












PW (Abinit) 


modi 6 


3.37 


14.8 


3.44 


13.9 


NAO (Siesta) 


modi'' 


3.38 


14.9 


3.44 


13.9 


OFDFT 












mcGGA 


spdl c 


3.37 


16.2 


3.43 


15.4 


TF+vW 


spdr 


3.37 


16.0 


3.43 


15.2 


GGA 


spdr 


3.37 


11.8 


3.46 


11.8 


TF+l/9vW 


spdl c 


3.37 


11.4 


3.46 


11.4 


OFDFT 












mcGGA 


modl d 


3.36 


16.2 


3.43 


15.2 


TF+vW 


modl d 


3.37 


15.9 


3.43 


14.9 


GGA 


modl d 


3.42 


10.8 


3.49 


10.1 


TF+1/9VW 


modl rf 


3.42 


10.3 


3.49 


9.5 



"Real space potential defined by Eq. (1201 ) (see text for details). 

''Real space potential defined by Eq. i'2\l ) (see text for details). 

^Reciprocal space potential defined by Fourier-Bessel transform 
of local potential Eq. 1201) and multiplied by f(q) function (see 
text for details). 

^Reciprocal space potential defined by Eq. (1 22 I t multiplied by f(q) 
function (see text for details). 



B. Pseudo-potential OF-DFT Tests 

1. OF-DFT Comparison for bcc Li 

For the OF-DFT bcc Li studies, we used a 128-atom 
supercell in Profess with the w sp di, Vmodi LPPs just 
described and both -E X c.lda and -E X c,gga- We did 
the Profess calculations for the Ttf, T TFvW ^=1,1/9, 
7gga,pbe-tw, and T mcG GA,PBE2 functionals. The com- 
puted E to t I atom values are plotted as a function of bcc 
lattice constant in Fig. [S] 

One sees that, as might be expected, the pure TF+XC 
model fails to bind. The pairing of other functionals, 
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Lattice constant (angstom) 

FIG. 8: Energy per atom vs. lattice constant for bulk 
bcc Li. OF-DFT results for Ttf, Ttf v w,o, Tgga,pbe-tw, 
and T mc GGA,PBE2 compared to the KS values. OF-DFT 
calculations with 128-atom supercell, VGGA.modi LPP, and 
-Exc,gga,pbe- KS calculations 2-atom unit cell with non- 
local PAW PBE pseudo-potentials (Vasp) and with Troullier- 
Martin PPs with PBE exchange-correlation, 2s 2 2p 2 basis set 
(8 NAO per atom) (Siesta). 

which we have discussed already, reappears. 7tfvW,c(=i 
pairs with T mc GGA,PBE2, and T TFvW)Q=1/9 pairs with 
Tqga.pbe-tw- The former pair gives a better descrip- 
tion of both the lattice constant and bulk modulus than 
the latter pair. The computed equilibrium lattice con- 
stants and bulk moduli are shown in Table UVl 

The equilibrium lattice constants predicted by the OF- 
DFT calculations with v sp di LPPs agree well with the KS 
PAW results. When the u mo di model pseudo-potential 
is used, the lattice constant from the OF-DFT calcula- 
tions with Tqga pbe-tw and T TFvW a=11 / g is an over- 
estimate of about 1 % for both LDA and GGA XC func- 
tional. This pair of OFKE functionals also predicts low 
bulk modulus values, again for both LDA and GGA XC 
cases. The mcGGA and TF+vW KE functionals do very 
well for the lattice parameter and slightly overestimate 
the bulk modulus value. 



2. OF-DFT Comparison for fee Al 

The utility of existing LPPs for OFDFT calculations 
obviously is a pertinent issue. To explore that, we con- 
sidered bulk Al. The model LPP in the form of Eq. 
(f!?2)) with parameters from Goodwin, Needs, and Heine 51 
was used in OF-DFT calculations. As before, this was 
done with the five OF-KE functionals, but here only in 
combination with the PBE GGA XC functional. Fig. [§] 
shows Profess results for a 4-atom fee cell compared 
with conventional KS results obtained with Vasp in the 
same cell with a 5 x 5 x 5 k-mesh calculation. The 




Lattice constant (angstom) 

FIG. 9: As in Fig. [8] for bulk Al but with a 4-atom super- 
cell. Siesta calculations performed with standard DZP basis 
set. The Goodwin, Needs, and Heine^ local model pseudo- 
potential used in the orbital-free calculations. See text. 

GGA and mcGGA KE functionals introduce numerical 
instability at expanded geometry. Aside from that, one 
again observes the same pairing of KE functionals as be- 
fore. The Tqga,pbe-tw an d Jt FvWiQ =i/9 functionals 
do not produce detectable minima. The X' mc GGA,PBE2 
and TrFvW,a=i pair predict equilibrium lattice constants 
(a = 4.05 and 4.06 A correspondingly), very close to the 
KS results (a = 4.05 and 4.09 A for PAW Vasp and 
NAO DZP Siesta calculations respectively). However, 
the shape of the two OF-DFT energy curves differs per- 
ceptibly from the KS results. In particular, the OF-DFT 
functionals predict a softer solid. 



V. SUMMARY DISCUSSION 

Several clear results emerge from this study. First, use 
of standard KS codes to solve the OF-DFT Euler equa- 
tion as a modified KS eigenvalue problem is problematic 
at best. At least for the all-electron case, it seems im- 
plausible as a productive route to routine OF-DFT calcu- 
lations. One could speculate that a better-behaved one- 
point approximate OFKE functional than mcGGA might 
not be such a challenge to standard KS algorithms. The 
repulsive nature of even the exact vg (recall Fig. [3]) makes 
that outcome seem rather doubtful. 

Second, even if a particular approximate one-point 
OFKE functional has singular behavior, it is possible 
that such a functional can deliver physically realistic re- 
sults. Those results can be obtained with a sufficiently 
refined direct Euler-Lagrange solution of the effective KS 
equation, Eq. ([2]). Thus, we are able to extract use- 
ful, self-consistent solutions for the recently developed 
simple mcGGA OFKE functional as well as the Tran- 
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Wesolowski GGA. These solutions enable understanding 
of the consequence of the singular behavior of their re- 
spective Pauli potentials. The Tran-Wesolowsk GGA has 
attractive singularities which cause strong over-estimates 
of the self-consistent density near the nuclear sites. In 
contrast, the properly positive mcGGA OFKE Pauli po- 
tential has positive singularities near the nuclei and the 
density is underestimated there. 

Third, we have presented a procedure for developing 
a local pseudo-potential for OFDFT calculations by do- 
ing a multi-channel weighting of a corresponding non- 
local pseudo-potential. The weighting is determined by 
KS calculations with the LPP such that the equilibrium 
non-LPP lattice parameter is reproduced. We showed 
that this yields a very good LPP. A remaining challenge 
for the OFDFT agenda is to construct a good LPP from 
an existing non-LPP without appeal to any bulk or ag- 



gregate system KS calculations. 

Fourth, once a suitable local pseudo-potential proce- 
dure is defined, the progress made on computational so- 
lution of the minimization problem for two-point OFKE 
approximations can be appropriated directly for use with 
one-point OFKE approximations. Even so, we do ob- 
serve numerical instabilities in the case of the mcGGA 
and GGA OFKE functionals. 
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